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ABSTRACT 

A propagation model of galactic cosmic protons through the Heliosphere was 
implemented using a 2-D Monte Carlo approach to determine the differential 
intensities of protons during the solar cycle 23. The model includes the effects 
due to the variation of solar activity during the propagation of cosmic rays from 
the boundary of the heliopause down to Earth's position. Drift effects are also 
accounted for. The simulated spectra were found in agreement with those ob- 
tained with experimental observations carried out by BESS, AMS and PAMELA 
collaborations. In addition, the modulated spectrum determined with the present 
code for the year 1995 exhibits the latitudinal gradient and equatorial southward 
offset minimum found by Ulysses fast scan in 1995. 
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Introduction 



During the last two decades - using balloon flights and space-borne missions -, the fluxes 
of Galactic Cosmic Rays (GCR) and their energy distributions were observed in different 
phases of solar activity. These data allow one to attempt a better understanding of pro- 
cesses related to the transport of OCRs through the Heliosphere. Furthermore, the study 
of propagation properties - i.e., the effect of solar modulation on the fluxes - of GCRs may, 
in turn, provide a tool to determine demodulated Local Interstellar Spectra (LIS) of GCR 
components, for instance, protons, light-nuclei, electrons, positrons, anti-protons, etc., thus. 
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addition, an accurate determination of demodulated sp ectra may allow one t o untangle fea- 
tures due to new physics - i.e., da r k matter (e.g ., see 
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. These spectra were measured i 



with an accuracy down to about or less than 30% and ii) covering a time duration longer than 
a solar cycle, i.e., these spectra were measured under solar conditions largely different. These 
data can be hopefully exploited to determine a general treatment of solar modulation in the 
inner heliosphere to be used for different phases of solar a ctivity and a better und erstanding 
of space radiation environment close to Earth (e.g., see iLeroy fc Rancoital 120071 . and refe- 



rences therein). In the near future even more accurate and systematic data will be available 
from AMS-02. This spectrometer is operational onboard of the Internation al Space Station 



from May 2011 an d is expected to collect data for more than a solar cycle ( iBattistoru 12010 



Bobik et al.ll2010al ). These observations will allow one to obtain accurate spectra with dif- 
ferent solar activity conditions from some hundreds MeV up to very high energy (a few 
TeV's); in addition, using the same experimental apparatus, systematic errors on measured 
fluxes are expected to be m inimized. Furthermore, observations made by the Ulysses space- 
craft ( ISimpson et al.lll992[ ) in the inner heliosphere could determine a latitudinal dependence 
of GCR (mostly prot ons) intensity with an equatorial so uthward offset minimum and a North 
polar excess (e.g., see lSimpson. Zhang and Bamelll996l ). Finally, it has be remarked that mo- 
dulation phenomena were observed at low ener gies (i.e., lower than 500 MeV/nucleon) in the 
ou ter heliosphere (e.g., s ee IWebber et al.ll2008l) and a re currently investig a ted, for in s tance , 
bv lLangner et all fcoosh : iLangner fc Potrieteil J2004h : iBobik et al.l J2008bl ): IPotrieteil J2OO8I ) 
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(see also references therein). 



In the present model, a two d imensional (2-D) - i.e., dependi ng on the helio-colatitude 
and radial distance from the Sun (IBobik et al.l |2003| . |2008| . l2010al ) - Monte Carlo approach 
is adopted to solve the transport equation of propagation of GCRs down to the inner he- 
liosphere, without addressing CR modulation observed in the outer. The model exhibits a 
slow time dependence because of the (almost) monthly averages of solar activity parame- 
ters adopted for the i) solar wind speed (Kw), h) tilt angle (at) of the neutral sheet and 
iii) diffusion parameter Kq (discussed in Sect. 12. ip . Furthermore, one has to remark that 
the solar wind usually takes of the order of or more than one year to reach the border of 
the heliosphere. As a consequence, the above parameters are locally evaluated within the 
heliosphere, allowing the modulation treatment to better (or dynamically) account for the 
effects of solar activity as a function of the distance from the Sun. In addition, the current 
treatment accounts for effects due to the charge sign of particles (i.e., the so-called parti- 
cle drift effect), e.g., those related, for instance, to a) the curvature and gradient of the 
interplanetary magnetic field (IMF) and b) the extension of the neutral current sheet inside 
the heliospher e. Thus the mod e l introduces a dep endence on the sign solar-field polarity 



{A) (e.g., see I Clem et al.l 1 19961 : iBoella et al.l 120011 ). The present code allows the fluxes of 



protons (as well antiprotons) and helium nuclei to be modulat ed from the border of the he- 
liosphere down to Earth - but outside Earth's magnetosphere (IBobik et al.ll2006l ) - in order 
to compare them with the available experimental observations. Furthermore, electrons and 
positrons modulated spectra can be deri ved accounting for the additional collision, radiative 
and inverse Compton energy- losses (see lBobik et al.ll201ld ). 



In the next sections, the heliosphere, drifts, diffusion tensor, determination of the dif- 
fusion parameter, dependence of both the solar wind and IMF on the radial distance and 
helio-colatitude, neutral current sheet are discussed (Sects. [2H1]). Then, the implementation 
of the mathematical model and the parametrization with the dynamical treatment of he- 
liosphere are treated (Sects. [5], [6]). Finally, comparisons among obtained modulated spectra 
of differential intensities with those experimentally observed are performed and discussed 
(Sects. EHU. 



2. Heliosphere and Drift Mechanisms 



The t ransp ort of galactic protons (CP) inside the heliosphere was initially treated 
by iParkerl (119651 ) . who demonstrated that - in the framework of statistical physics - the ran- 
dom walk of the cosmic ray particles is a Markoff process, descr i bable by a Fokker-Pl anck 
equation (hereafter FPE) (e.g., see also lAxfordlll965l : lFisklll976l : iPotgieter et al.lll993l . and 
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also Sections 4.1.2.4 of iLeroy fc Rancoita 2011 . and references therein). Thus (at the time 
t), the number densit}0 U of GPs per unit interval of particle kinetic energy T (the so-called 
differential density) can be obtained from the solution of the FPE: 



dU 
'dt 



d_ 
dxi 



^ = ^ K. 



dU 



1 dV. 
+ 3 



d 



dxi ''^'^ 



(e.g., see IJokipii et al.l 119771 . Equation (4.75) in Section 4.1.2.6 of iLeroy fc Rancoital 12011 
and references therein) with V^^.i the solar wind velocity along the axis Xi, 



Vd,, 



dXj 



(2) 



the dr ift velocity (e.g., see IJokipii et al.l 119771 : IJokipii fc Levyl 119771 and also iBobik et al. 
2010bl and references therein), 
diffusion tensor - respectively - 



2010bl and references therein), and K^, the antisymmetric and symmetric part of the 



T + 2m^c2 



T + rrirC^ 

and rrir the rest mass of the proton. The number density U is related to the differential 
intensity J as: 



where v is the speed of the GCR particle. Equation ([T]) - as well known - describes i) 
the diffusion of GCRs by magnetic irregularities, ii) the so-called adiabatic- energy changes 
associated with expansions and compressions of cosmic radiation, iii) the convection effect 
resulting from the solar wind with velocity and iv) the drift effects related to the drift 
velocity (vd)- In turn, the drift velocity is determined by the antisymmetric part of the 
diffusion tensor [see Eq. (jS]) and Sect. H] which accounts for gradient, curvature and current 
sheet drifts of particles in the IMF, i.e., it depends on the charge sign of particles. 



Furthermore - as discussed by IJokipii &: Levyl (119771 ) -, one can re- write Eq.([T]) as 



dU 
'dt 



d_ 

dxi 



dXj 



1 dV^^i d 



3 dx, dT 



(a^eiT U) 



d_ 

dxi 



[(Kw,i + Vd,i) U] . 



(4) 



Thus, one obtains that drift effects are accounted for by a convection velocity in which the 
drift velocity is added to the solar wind velocity. In this way, the resulting effective convection 



^The equivalent expression in terms of the omnidirectional distribution function of CR part icles with 
mome ntum p, at the position rand time t can be found, for instance, expressed in Equation (1) of lPotgietei 
(|l998t ) (see also references therein). 
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velocity may n o n-neK ligibly differ from tliat due to tlie solar wind; but - as remarked by 
Jokipii &: Levyl ( 119771 ) - noting tliat V ■ = 0, one finds that drift effects do not contribute 
to tlie adiabatic-energy clianges [second riglit-fiand term of Eqs. ([H H])]. Even i f drift effects 



are 



i ncludeci§ in Eqs. (HI H]), some modulation model^ neglected it (e.g., see Jokipii et al. 



19771 : lUsoskin et al.ll2005l . and references therein). Gradients of particle density can also result 
from the convection effect. Drift mechanisms can modify both the radial and (solar) latitude 
dependence of the gradient magnitude. For instance, drift motions can affect modulated 
GCR spectra by redirecting particles within the heliosphere (jjokipii et al.lll977l ). When the 
particle Larmor radius is much shorter than the magnetic-field scale length, drift effects can 
be taken into account by evaluating the average distance in which a relevant field variation 
occurs. Drift effects affect particle motions over large distances due to the large scale variation 
of the IMF strength. Different intensities of GCR modulation were ob s erved in time periods 



with o p posite field polarit y, for instance , bylEmerson &: Meyerl ( ll984j ): iGarcia-Munoz et al. 



( 119861 ): IClem et al.l ( l2000l ): iBoella et al.l (l200l[ ). Thus, it is necessary to explicitly consider 
particle drifts inside the equation of propagation of GCR. 

As well known for a reference system with the 3rd coordinate along the average magnetic 
field, the symmetric part of the diffusion tensor (or coefficient) - for an isotropic perpe ndicular 
diffusion - includes both the transverse (K \ ) and parallel (K \i) components (e.g., see I Jokipii 
197ll : iPotgieter fc Moraallll985t iPotgieter &: Le Rouxlll994j ). In turn, for a standard Parker 
field [Eq. ( 1T5|) ] these two components are related to the radial component in heliocentric 
spherical coordinates as 

Krr = Kucos'^ip + K±sin'^ip, (5) 



with ip the angle be tween radia 



[Eq. ([IS])] - (e.g.. seelFisk 



1976 



and magnetic field direct ions - the so-called spiral angle 
Potgieter &: Le Rouxlll994l ) and Kgo = K^, where 6 is the 



polar angle ( Potgieter et al. 1993[ ). It has to be remarked that the general transformations 



of the symmetric and antisymmetric parts of the diffusion tensor fro m field-aligned to helio- 
spheric (spherical) coordinates can be found in ( iBurger et al.ll2008l ). Furthermore, it has to 
be remarked t hat a general discussion a bout the role of parallel and perpendicular diffusion 



is available in iGiacalone fc Jokipiil Il999 



Potgieter fc Le Rouxl ( 1l994l ) (see also IPotgieter et al.l Il993l ) suggested that the parallel 



Parker 


1965; 


Jokioii & Parker 


1970l 


Jokipii & Lew 


1977; 


Jokipii et al. 



19771: |Potgieterj|l998t) . 



•^Like, for instance, the so-called force-field model (FFM) (see Gleeson &: Axford 19681) . 
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diffusion coefficient is given by 



K\\^ (3h{r,t)Kp{P,t) 



3B 



(6) 



with (3 = v/c , V the particle velocity and c the speed of hght; the diffusion parameter ki 
accounts for the dependence on the solar activity and is treated in Sect. I2.1j (typically 
pa 5 nT) is the value of IMF at Earth's orbit but it varies as a function of the time; B is the 
magnitude of the large scale IMF (discussed in Sect. |3]), thus, it depends on the heliospheric 
region (Sect. [5]) through which GCRs are transported; finally, the term Kp takes into account 
the dependence on the rigidity P of the GCR particle and is usually expressed in GV. To a 
first approximation, one can assume that 



for particle rigidities above a threshold value Pth. 



commonly suppos e d by many authors (e. g., see 



19681; IPerkol Il987l : IPotgieter fc Le Rouxl Il994 



(7) 



within the rigidity range (0.4-1.015) GV, as 



Gloeckler &: Jokipiilll966l : iGleeson &: Axford 



Strauss et all 1201 ll ). In the present model. 



Kp is assumed to be equal to the value of the rigidity (P) above the upper limit of the 
th range, i.e., for proton kinetic energies > 0.444 GeV (see Sects. 



Pth range, i.e., for proton kinetic energies > 0.444 GeV (see Sects. 17.21 17.2. ID. Below Pth, 
it ca n be usually approxima t ed to a constant (e .g., see iPerkd Il987l : IPotgieter fc Le Roux 



I994J : IWibberenz et al.l l2001t IStrauss et al.l I2OIII ) . It has to be remarked that nowadays 
treatments resulting in a more comple x dependence of th e diffusion tensor on rigidity are 



proposed by several authors (e.g., see iFerreira et al.l I2OOII ; iPei et al.l l2010al and references 



therein). Some of these studies are motivated from dealing with magnetohydrodynamic tur- 
bulence in the expanding solar wind and/or accounting observations carried out on data 
of low energy electrons collected u sing spacecrafts [for instance, (3-1 0) MeV from Ulysses 
space craft in (IFerreira et al.l 1200 ll ) and 16 MeV from Pioneer 10 in ( IPotgieter fc Ferreira 



2OO2I)] 



In heliocentric spherical coordinates, the perpendicular diffusion coefficient has two com- 
ponents, one along the radial direction, K_i_r, the other one for the polar direction K±g. 
is the ratio between perpendicular (in the radial direction) and parallel diffusion coeffi- 
cients, i.e., K±r = PkK\\. In t he present mo del, we u se Pk = 0.05: this value is in the 



mid of the range s uggested by iPalmerl (119821 ) (see also iGiacalond Il998l and Section 6.3 of 
Burger et al.ll2000l ). The value of the perpendicular diffusion coefficient in the p olar direc 



tion (K I 0) can be assumed to be aln aost equal to th at radial {K±r) (e.g., see IPotgieter 
2OOOI . and references therein). However, IPotgieter J2OO0I ) suggested the usage of an enhanced 
K±9 in the polar regions in order to reproduce the amplitude and rigidity dependence of 
the latitudinal gradients of GCR differential intensities for protons and electrons (e.g., see 
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Potrieter 


1997; 


Heber et al. 


1998) 



tion, e.g., see Figure 7 in that article) in the colatitude regions 120° < 6' < 130° and 
60° > 6 > 50°. He also derived that K i g has to be increased by a factor of about (or larger 
than) 10; iFerreira fc Potgieter used a factor of 8. In the current code, Kj_g is given 

by: 



K 



±9 



10 K±r, in the polar regions, 
Kj_r, in the equatorial region. 



where the polar regions correspond to colatitudes with 6 < 30° or 6 > 150°, while the 
equatorial region corresponds to colatitudes with 30° ^ < 150°. The solar colatitudes of 
30° and 150° correspond to those at which the SW speed becomes constant in periods not 
dominated by high solar activity [Eq. (1171) ]. The usage of the transition function can be 
fully implemented in current treatment, but is not required with the present overall code 
accuracy. In fact, the results obtained from the so-called "L" model (i.e., the one with a better 
agreement with data, see Sects. 17.21 17.2. ip indicate that only in periods not dominated by 
high solar activity the enhancement of K^g [Eq. ([8])] slightly improve the overall agreement 
with data by a few percent. Finally, in Appendix |A] the diffusion coefficients in heliocentric 
polar coordinates are expressed in terms of those parallel and perpendicular to the IMF. 

Moreover, it has to be remarked that the diffusion tensor i) is not well determined during 
solar maxima and ii) can be adapted to better account for the complex structu re of the IMF 



- whi ch depends on the solar activity - found with Ulysses spacecraft (e.g., see iBurger et al. 
20081 . and references therein). For instance, Potgieter, Burger and Ferreira (2001) - see also 
references therein - discussed the so-called propagating diffusion barriers and suggested a time 
dependent model for the diffusion coefficients. The latter are supposed to be oc [BQ/B(t)]"-, 
where B{t) is the IMF magnitude at the time t an d Bo = 5 nT is the average IMF magnitude 



durin g minimum modulation conditions at Earth (jPotgieter fc Ferreirall200ll ; iPotgieter et al. 



20031 ) and n is the ratio betw 'een the actual tilt-angle value (Sect. [3l) and that close to solar 



minimum (7°-15°) (e.g., see IPotgieter fc Ferreira! l2001t IPotgieter et al.l 120011 ). However, in 



the current model the time dependence of diffusion coefficients is taken into account using a 
diffusion parameter, which is treated in Sect. 12.11 The agreement with data obtained during 
high solar activity is discussed in Sect. 17.21 



2.1. Diffusion Parameter in the Framework of the Force Field Model 



In the F FM (e.g.. see iGleeson fc AxfordI Il968l ; iGleeson fc UrchI Il97ll and also Sec- 
tion 4.1.2.4 of lLeroy &: Rancoitall201ll ). Gleeson & Axford (1968) assumed that, at the time 
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t, i) modulation effects can be expressed with a spherically symmetric modulated differen- 
tial number density U of GCRs, ii) the diffusion coefficient reduces to a scalai0 given by a 
separable function of r (the radial distance from the Sun) and P (the particle rigidity in 
GV): 

IC{r,t) = l3h{r,t)Kp{P,t) (9) 

with Kp from Eq. ([7]) for particle rigidities above ~ 1 GV, and iii) the modulation occurs in 
a steady-state condition, i.e., the relaxation time of the distribution is short with respect to 
the solar cycle duration so that one can assume that the partial derivative of U with respect 
to time is zero. They derived that the differential intensity [Eq. ([3])] at a radial distance r is 
given by the expression 



J{r,Et,t) = J(rto,Et + $F 



2^4 



[(Et + $p)2-m: 



2„4 



(10) 



where J(rtm, Et + $p) is the undisturbed intensity beyond the solar wind termination lo- 
cated at a radial distance rtm from the Sun; is the total ene rgy of the particle with 
rest mass rrir and, fina lly, $p is the so-called force-field energy loss f Gleeson fc Axford 19681: 
Gleeson fc Urclj|l971l) . When mod ulation is small [i.e., for $p ^ rrij-c^, T] ( IGleeson fc Axford 



1968 



Gleeson fc UrchI Il97ll . Il973[ ) , they determined that 



ZeP 
KpiP,t) 



(^s{r,t) ^ Ze(f)^{r,t), 



where Ze is the particle charge and (f)s{r,t) is the so-called modulation strength (or modu- 
lation parameter) usually expressed in units of GV (or MV). Assuming that Vsw (the solar 
wind speed) and k^ are almost constant, .s(^, ^) is linearly dependent on (rtm — r) (e.g., see 



Equation (4.64) of lLeroy fc Rancoitall201ll ). from which one gets that the diffusion parameter 
is given by 

^^^'^ 3Mr,t) ' ^^^^ 

i.e., ki [similarly to (f)s{r,t)] is linearly dependent on (rtm ~" ''")• As already mentioned, in the 
FFM the diffusion coefficient /C(r, t) is a scalar quantity and does not account for effects 
related to the charge sign of the transported particl es. (l)Jr,t) is independen t of the species 
of GCR particles (e.g., see discussion at pag e 1014 oflGleeson &: Axfordlll968l or Equation (1) 



of Usoskin and collaborators 2005, and also Bobik et al 



2011a 



3). Usoskin and collaborators 
2005 monthly determined the values of the modulation strengths [0s(^Earth)] ioi the time 



"'While in Eq. , it is expressed by a tensor with a symmetric and an antisymmetric part (see discussion 
in Sect.[2|). 
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period from 1951 up to 2004 using measurements of neutron monitors (i.e., located at rEarth = 
1 AU); while the values of solar wind speeds are available from NASA/GSFC's OMNI data 
set through OMNIWeb. 
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Fig. 1. — Diffusion parameter Kq (left side) and percentage differences -Rperc [Eq. ( |T^ ] (right 
side) as a function of the SSN value; the central continuous lines are obtained from a fit of 
Kq with respect to SSN values in the range 10 < SSN < 165; the dashed and dotted lines 
are obtained adding (top) or subtracting (bottom) one standard deviation from the fitted 
values. 

To determine 0s('^Earth)5 Usoskin and collaborators (2005) used an approximated ex- 
pression of the Local Interstellar Spectrum (LIS) for protons from Burger, Potgieter and 
Heber (2000). In practice, their spectrum differs from that due to Burger, Potgieter and 
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Solar polarity 


A < 


A < 


A > 


A > 


and phase 


ascending 


descending 


ascending 


descending 


Cl 


+0.0001686 


+8.872x10-5 


+2.39708x10"^ 


+2.28037x10-"^ 


C2 


+0.001488 


+0.001874 






C3 






-8.28987x10-^ 


-1.00984x10-6 


C4 


-3.164x10"^ 









Table 1: Parameters of the polynomial expression ( IT3|l . 



Heber (2000) by about or more than 5% at kinetic energies lower than about 117 MeV. Fur- 
thermore, in the present work we found that the error-weighted average of the differential 
spectral index, 7„a [Eq. (!3T]) and discussion in Sect. 17.1] , of the proton LIS is only compa- 
tible, within one standard deviation, with the differential spectral index (7 = 2.78) of the 
spectrum from Burger and Potgieter (1989) or Burger, Potgieter and Heber (2000). It has 
to be remarked that the latter spectral index is the one used by Usoskin and collaborators 
(2005). Usoskin and collaborators (2005) [see Appendix A in that article] also found that 
using other commonly adopted LIS's their corresponding values of the modulation strengths 
follow a linear relation with respect to 0s('"Earth)- However, the differential spectral indexes of 
these spectra are not compatible within three or more standard deviations with that found 
in Sect. 17.11 Moreover, it has to be noted that the response of neutron moni tors has to be 



evalu ated by combining a) the effects of both the geomagnetic cutoff rigidity (lUsoskin et al. 



20051 ) which results in a r educed sensitivity of d etection apparata and b) the so-called at- 



mospheric yield function (jClem fc Dormarul2000l ). Thus, one finds that i) the contribution 
of the GCRs with rigidities below 2GV amounts to about or less than 1.1% of the to- 
tal neutron monitor counts due to particles with energies up to about 50 GV and ii) the 
maximum of neutron monito r sens itivity - i.e., the maximum of the response function [see 
Figure 7 of I Clem &: DormanI (120001 )] - occurs in the rigidity interval (3-15) GV. In addition, 
Boella and collaborators (2001) determined - using IMP8 satellite data during the period 
1973-1995 - that charge effects (discussed in Sects. [H [2]) result in a variation of proton or 
helium fluxes during solar minima with opposite magnetic field polarities of 14 ± 6% at 
300 MeV/nucleon. Th is variation steadily decr eases with increasing energy [e.g., see Fi- 
gure 4.13 at page 378 of iLeroy fc Rancoital (I2OIII )]. As a consequence, 0s('"Earth) is expected 
to be marginally affected by drift effects. 



ki [Eq. ( ITT]) ] depends on the value of solar wind termination locat ed at a radial distanc e 
related, in turn, also to solar wind speed [e.g., see Chapter 7 of iMeyer- Vernetl ( 120071 ). 
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99 AU 

h- = 99 AU 

(^tm ~ ''"Earth ) 



(12) 



and Sections 4.1.2.3, 4.1.2.4 of iLeroy fc Rancoital ( 120111 )]. In the present simulation code, 
the effective heliosphere assumes that the solar wind termination is located at 100 AU (see 
a further discussion in Sects. I7.3p . Therefore, from the diffusion parameter ki one has to 
derive that (Kq) for an effective heliosphere with a radial extension of 100 AU. In practice, 
for a radial extension of 100 AU the diffusion parameter Kq [Eq. ( IT2] replaces ki in Eq. ([6]) 
(for instance, see Appendix and in allows one to obtain similar modulation effects on the 
differential intensities of GCRs with respect to those obtained using ki when the heliosphere 
has a variable radial extension rt^. Using Eq. ( 1111) one obtains 

.3 0s(rEarth) 

where 99 AU (as already mentioned) is the distance of the Earth from the border of the 
effective heliosphere used in the current simulation code. In Fig. [H the diffusion parameter 
Kq - obtained from Eq. (fl2|) - is sho wn as a function of the corresponding value of smoothed 
sunspot number, SSN, (ISSNI I2OIOI ). The Kq data had to be subdivided in four sets, i.e., 
ascending and descending phases for both negative and positive solar magnetic-field polari- 
ties. For each set, the data could be fitted with a practical relationship (see Fig. [1]) between 
Ko and SSN values for 10 < SSN < 165, i.e., finding 

Kf = ci + C2X SSN"^ + C3 X SSN + C4 x SSN^ (13) 

with the parameters Cj shown in Table [H In addition, the data were found to exhibit a 
Gaussian distribution of percentage differences (-Rpcrc) of Kq values from the corresponding 
fitted values Kp, with 

The rms values of the Gaussian distributions were found to be ~ 0.1339, 0.1254, 0.1040, 0.1213 
for the phases ascending with A < 0, descending with A < 0, ascending with A > 0, descend- 
ing with A > 0, respectively. From the practical relationship found [Eq. ( 1T3|) ]. we can use 
the estimated SSN values to obtain the diffusion parameter Kq at times beyond 2004. This 
procedure allows one to extend the ~ 40 years period by exploiting the practical relation- 
ship between the fitted Kq values and the SSN values (one of the main parameters related 
to the solar activity). In addition, we introduced in our code a Gaussian random variation of 
Kq with rms's corresponding to those found for each subset of data. Results of the simula- 
tion with and without the Gaussian variation are consistent within the uncertainties of the 
code. Furthermore, it can be noted that Kq results in providing an overall increasing (for 
rtm lower than 100 AU) or decreasing (for rtm larger than 100 AU) of modulation effects. A 
tuning of the effective extension of the heliosphere and its dependence on the solar activity is 
likely to be obtained using the experimental data from long-duration accurate observations, 
like those from the AMS-02 spectrometer. 
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3. Solar Wind and Latitudinal Dependence IMF 



Parker (1958) suggested that the solar corona is stationary expanding due to an outflow 
of the coronal plasma - generating the so-called solar wind - with a spherically symmetric 
velocity. In his model, the solar wind speed becomes almost co nstant (Kw) b eyond a radial 
distance from the Sun r{, ^ (0.3-0.4) AU (e.g., see Figure 1 of |Parkerlll958l ). Furthermore, 
the magnetic-fleld lines are frozen in the streaming particles of which the solar wind con- 
sists. Thus, beyond r^, in a spherical reference frame rotating with the Sun the components 
of the outward velocity of a plasma element carrying the magnetic fleld are: K = V^w, Ve = 
and = uj{r — rf,) sm6 with u the angular velocity of the Sun. The streamline has the shape 
of an Archimedean spiral (termed Parker spiral). 



In heliocentric spherica l coordinates, the standar d Parker spiral fleld can be expressed 
as (e.g., see Equation (2) of iHattingh &: Burgerlll995l ): 

A 



{cr-Ve^) [l-2H{e-e')] 



(15) 



where y4 is a coefficient that determines the fleld polarity and allows \Bp\ to be equal to 
(Sect. [2]), i.e., the value o f IMF at Earth's orbit as e xtracted from NASA/GSFC's OMNI data 
set through OMNIWeb (IKing fc Papatashilill2005l ): and e^j, are unit vector components in 



the radial and azimuthal direction, respectively ; 6 is the co-latitude (po lar angle); 6' is the 
polar angle determining the position of the HCS (jJokipii fc Thomaslll98ll ); H is the Heaviside 
function, thus, [1 — 2H{6 — 6')] allows Bp to interchange the sign in the two regions - above 
and below the heliospheric current sheet (HCS) - of the heliosphere; finally, 

u (r — Tft) sin 6* 



r = tan ip 



(16) 



with ip the spiral angle. In the present, model u is assumed to be independent of the helio- 
graphic latitude and equal to the sidereal rotation at the Sun equator. However, the simple 
representation of the Parker spiral [Eqs. f fT5| [T6l) ] based on a constant solar wind speed 
needs to be complemented with the present knowledge of the speed [l^sw(6')] dependence on 
solar colatitude. Large variations of the solar wind structu re were observed for solar lati- 
tudes up to 1 80° I by Ulysses spacecraft (IWenzel et al.lll992l ). For instance, during a period 
of low solar activity the solar wind speed increases by almost a factor two from the eclip- 
tic pla ne to poles, thus subd ividing the heliosphere in two regions with slow and fast solar 
wind (IMc Comas et al.l 120001 ). For representing the observed speeds, Fichtner, Ranga and 
Fahr (1996) suggested that the solar wind speed may be proportional to (1 + cos^ 9). In the 
present model we use: 



K^'.nax, for 6 < 30° and 6 > 150° 



;i + Icos^l), for 30° <9 < 150° 



(17) 
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with Vswn,!,^ — 760km/s (e.g., see lMc Comas et al.ll2000l )] and Vsw^i,-, is t he corresponding va- 
lue ex tracted from NASA/GSFC's OMNI data set through OMNIWeb Jxing fc Papatashih 
20051 ) . Equation f[T7|) exhibits a shghtly better agreement with observed data than that pro- 
posed by Fichtner, Ranga and Fahr (1996). Jokipii and Kota (1995) and Pommois, Zimbardo 
and Veltri (2001) proposed other functions for such periods. However, these functions de- 
pends on an additional parameter related to the latitudinal extension of the region with a 
slow solar wind. The parameter can be determined only using measurements to be performed 
largely outside the ecliptic plane, hke those due to Ulysses spacecraft. Thus, Eq. (IT7|) has the 
advantage to allows one to more generally treat periods of low solar activity. Furthermore, 
McComas and collaborators (2000) observed that during the Sun's approach to solar maxi- 
mum a) the coronal structure becomes increasingly complex and b) the magnetic field be- 
comes less dipolar. In the present model, for the solar wind we assume a speed independent of 
the colatitude in periods characterized by a large solar activity. As previo usly, the speed value 
is ex tracted from NASA/GSFC's OMNI data set through OMNIWeb Jxing fc Papatashih 
2005h . 



Potgieter et al.l ( 1l989l ) pointed out how classical drift modulation models - based on 
the Parker magnetic-field up to the polar region - encounter difficulties (see also Sect. 17.41) 
i n acc ounting for the significantly lower latitudinal dependence of CRs intensity. ISimpson 
(119961 ) subsequently observed this phenomenon using Ulysses spacecraft data collected in the 
inner heliosphere. Heber and collaborators (1998) remarked that a) one needs to assume an 
anisotropy of perpendicular diffusion coefficient and enhancement in the latitude direction 
(as already treated in Sect. [2]), and b) Parker's IMF has to be modifiecj^ as proposed by 
Jokipii and Kota (1989). 

In the present model, the magnitude of th e magnetic field [Eg. ([TSD] is enhanced int ro- 
ducing a small latitudinal component (e.g., see lLangneiil2004l : iLangner fc Pot gieted 120041 ) 



with Tq the solar radius. 



Bo 



A 



m 



8.7 X 10" 
sin 6 



(19) 



for e < 1.7° and 6 > 178.3°, 5 is ~ 3 x lO^^ flFichtner et al.lll996r ). It has to be noted that 
Eqs. flTHl [TI?]) allows one to obtain V ■ -B = 0. The magnitude of the magnetic field used in 



^Limited to polar regions, Fisk (1996) proposed that magnetic-field lines are non radially expanding. In 
addition, Hitge and Burger (2010) have attempted to merge Parker and Fisk magnetic-fields into a hybrid 
field. 
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Fig. 2. — Magnitude {\B\) of the IMF (dashed hne) from (IKota &: Jokipiilll989h - computed 
using Eqs. f lTB] - [2U]) - compared with that from Parker (sohd hne) [Eq. fll5p ] at 1, 5 and 
10 AU as a function of the colatitude. For the purpose of this calculation, at lAU and 90° 
\B\ = 5nT. 



the current model is given by (jJokipii &: K6talll989l ): 



B 




r2 + 



r 

re 



(20) 



In Fig. [21 the magnitude (|-B|) of the IMF from (iKota fc Jokipiil Il989h is compared with 
that from Parker [Eq. ( ITSl) ] at 1, 5 and 10 AU as a function of the colatitude: the field 
magnitude significantly increases in the polar regions (colatitude < 10° and > 170°), while 
it is almost unchanged in the ecliptic region (colatitude ~ 90°). \B\ (Fig. |2]) was computed 
using Eqs. f|T8] - [20|) . As discussed by Haasbroek and Potgieter (1995), the above modification 
of the Parker IMF allows the modulation effect in the polar regions to be increased and, 
subsequently, more reahstic radial and latitudinal gradients. 



As well known (e.g., see lBravo et al.lll998l and references therein), during several years 
around solar minimum the general structure of the solar magnetic-field is more or less axi- 
ally symmetric, dominated by the dipole component. These periods are char acterized with 
corresponding low values of the tilt angle (at < 30°, IPotgieter et al.ll200ll ). As the solar 
activity increases, the dipolar structure inclines more and more with respect to the ro tation 
axis and the effect of higher multipoles becomes more relevant ([Sanderson et al.l [2003[) . Du- 
ring the years of high activity, the structure of the solar magnetosphere is very complex 
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and the dipole component is very tilted (e.g., see ISanderson et al.l l2003l : IWang fc Sheeley 



2OO2I ). Th ese periods are char acterized with corresponding large values of the tilt angle 



at > 75°. |Potgieter et al.ll2001). Fin ally, one can remark that, dealing with neutron-monitor 



measurements, iCliver fc Lingl (120011 ) concluded that a diffusion/convection dominated mo- 
dulation occurs when the tilt a ngle exceeds 50° [for a previous discussion and illustration 
with numerical models, e.g., see iPotgieterl ( 1l995l ) and references therein]. 

In the current model the evolution of the solar magnetosphere and, subsequently, of the 
IMF is (partially) taken into account using the diffusion parameter (treated in Sect. 12. ip and 
the actual value of tilt angle. In fact, the diffusion parameter depends on the solar phase (as- 
cending and descending) and solar polarity (positive and negative), and is practically related 
to the actual value of smoothed sunspot number via Eq. f lT^ : while the tilt angle allows one 
to gradually modify the contribution of drift effects to modulation (Sect. H]). The agreement 
with data obtained during high (low) solar activity is discussed in Sect. 17.21 (Sect. 17.2. ip . 



4. The Neutral Sheet and Large Scale Gradients of IMF 



Ka expresses the value of the antisymmetric part of the diffusion tensor and results from 
the effects on the motion of cosmic-ray particles due to drift mechanis ms. In a coordinate 
syste r n with the 3rd coordinate along the average IMF, one finds (e.g., see lPotgieter fc Moraal 
19851 : burger k Hattinghlll995h 

Ka = TT^, (21) 
3 Ze\B\ 

where p, v and Ze are the momentum, velocity and charge of the cosmic-ray particle, respec- 
tively. Thus, the antisymmetric elements of the diffusion-tensor matrix (Sect. |2]) are 



with ejj^fc the Levi-Civita symbol (e.g., see Equation (10) of |Parkerlll965l ). 

As already mentioned in Sect. El the drift velocity va [Eq. (E])] accounts i) for effects due 
to gradient and curvature drifts experienced by cosmic-rays particles transported trough the 

■ y (e.g., see 



Parker 


1957; 


Burger. Moraal & Webb 


1985; 


Potgieter & Moraal 


1985) and iii) 


lated using the antisymmetric part of the diffusion tensor (e.g., see 


Parker 


19651 



19771 ; iPotgieter &: Moraallll985l ; iBurger &: Hattinghlll995l . and references therein) 



Burger, Potgieter and Heber (2000) [see also references therein and (jPalmerl Il982 



Lockwood fc Webbeiil2000l )] remarked that the observational results (carried out below 5 GV) 
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are consistent with a small (or very small) ratio of the perpendicular to parallel diffusion 
coefficients. As discussed by Parker (1965), a small value of that ratio indicates that cosmic- 
rays particles are practically moving through several gyro-orbits between each scattering 
event, i.e., drift motion is weakly affected by scattering. In addition, for cosmic-rays particles 
with rigidities < (10-15) GV and an IMF expressed by Eqs. (fT5l[20!) . the particle gyro-radius 
is smaller (or much smaller) than any local (i.e., inside the heliosphere) scale variation of 
magnetic field L = \{l/B){dBi/dxi)\~^. In this way, for regions outside that of HCS, Isen- 
berg and Jokipii (1979) remarked that Vd is determined by the terms due to the gradient 
and curvature drifts (e.g., see also iParkerl 119571 : [Armstrong et al.lll985l ). 



Potgieter and Moraal (1985) treated the modulation of GCR's for steady state condi- 
tions with relevant drift effects including that due to a wavy HCS (WHCS). They succeeded 
in formulating a 2-D description (of the WHCS), which - as discussed by Burger and Hat- 
tingh (1995) - is equivalent to the treatment of transport in a three-dimensional heliosphere 
with the assumption of an axis-symmetric particle distribution. Thus, they allowed one to 
neglect the azimuthal dependence. The effect of a WHCS was included via an appropriate 
modification of the antisymmetric part of the diffusion tensor. In this 2-D modeling, the 
WHCS is described as a wide region whose width depends on the rigidity of cosmic-ray 
particles and actual value of the tilt angle (at)- The resulting drift velocity in heliocentric 
polar coordinates - as u sed in the current model - is given by (e.g., see Equation (6) of 
Burger fc Hattinghlll995l ): 



Vd 



fi9)V x{Ka^] + 



dm 

dd 



Ka 
r 



X 




(22) 

= ^^dr + ^^HCS, (23) 

where Ka is from Eq. fl2T]) . 6 is the colatitude, f{6) is a. trans ition function that accounts for 
the effects of a wavy neutral sheet (jPotgieter &: Moraallll985l ) and eg is t he unit vector along 
the la titudinal direction. /(6') is expressed as (e.g., see Equation (14) of iPotgieter fc Moraal 
1985h : 

{{l/ah) arctan{{l - [(2 6')/7r]} tan(a/,)} , if < |, 
l-2H[e-{'n/2)], ifc, = f 
with H the Heaviside function, 



arccos 



TT 



2ch 



- 1 



e.g., see Equation (15) of IPotgieter fc Moraallll985l ). 



TT 1 

Cfe = - - - sin(at + A6'hcs) 
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e.g., see Equation (23) of iBurger fc Potgieterlll989l ). 



2r, 



HCS 



(rp is the particle gyro-radius, e.g., see lHattingh fc Burgerlll995l and also Section 4.2 of lBurger fc Hattingh 

19951 ). finally, f{ch) = 0.5 and /(7r/2) = 0. A6'hcs is determined from the maximum distance 

that a particle drifting along the neutral sheet can be away from this sheet (IBurger &: Potgieter 

19891 ). The first term (-Udr) of Eq. ( I22l) accounts for the gradient and curvature drifts, 

the second ("yHcs) for drift in the region affected by a WHCS. The transition function 

sets t he rate at which the firs t term of Eq. fl22|) goes to on the ecliptic plane {6 = 

■k/2) JPotgieter k MoraallEgSsf ). 



5. Parameters of the Effective Heliosphere used in the Current Model 



As discussed by Potgieter (2008) (see also references therein), until recently the helio- 
sphere was assumed to be spherical in most modulation models with an outer boundary 
at radial distances beyond ~ 100 AU. Presently, the heliospheric structure is considered la- 
titudinally asymmetric (particularly) during solar minimum conditions mostly because the 
SW depends on the latitude and solar activity (Sect. [3]). As a consequence, the position of 
termination shock (where the SW ram pressure is balanced by interstellar pressure), TS, can 
exhibit a latitudinal asymmetry. 

Using solar wind speeds observed from Uly sses, Whang and collaborators (e.g., see 
Whang fc Burlagalbood : Iwhang et al.lboosl . 120041 ) could estimate the radial position of TS 
on and outside the ecliptic plane. They found that a) on the ecliptic the radial distance 
of TS is about of 80 AU on average (without large variation between low and high solar 
activities), b) near the ecliptic the radial distance varies by less then 20 AU and c) outside 
the ecliptic plan e (e.g., at a latitude of 35°) the location of the TS increases by more than 
or about 50 AU (IWhang et al.ll2003l ). In addition. Whang and collaborators estimated that 
the averaged value over a 26-years period o f the radial distance of the TS increases with 
latitude [see Table 2 of (IWhang et al.l l2003l )]. It is worthwhile to remark that ^ 100 AU 
is the averaged value ov er the correspondin g solid angle of the TS location, which can b e 
obtained from Table 2 of (IWhang et al.ll2003l ). Furthermore (e.g., see lStone et al.ll2005l . l2008l ). 
Voyager 1 and 2 reached the TS in 2004 and 2007 located at about 94.0 AU and 83.7 AU, 
respectively, in agreement with the predictions from Whang and collaborators. Langner and 
Potgieter (2005) treated symmetric and asymmetric TS models and concluded that for A > 
cycle for solar minimum no significant difference occurs; for A < cycle differences remain 
insignificant in nose direction while, approaching the tail direction, some differences can be 
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appreciated at proton energies below (1-1.5) GeV. However, Langner and Potgieter (2005) 
and Potgieter (2008) suggested that, in general, a symmetric TS with a radial distance of 
~ 100 AU is still a reasonable assumption. 



In the present model (as already discussed in Sect. 12. ip . the effect of the modulation is 
obtained for the GCRs propagation trough a symmetric effective heliosphere with a radius 
of 100 AU. The diffusion parameter Kq is determined (followin g the proc edure described in 
Sect. 12. ip using the values of modulation strength, SSN values ( 1SSN||2010| ) and radius of the 
effective heliosphere. Furthermore, it has to be remarked that (see discussion in Sect. 12. ip the 
atmospheric yield function results in a diffusion parameter related to modulated intensities 
of GCRs (mostly protons) with rigidities above 2 GV. 

Other parameters (which depend on the solar activity) are the tilt angle at of the 
HCS, magnetic field polarity [related to the sign of the coefficient A in Eq. f lT^ ]. mag- 
netic field magnitude (-B®) and solar wind velocity (V^w)- The latter two parameters are 
measured at Earth's orbit. The polarity of the magnetic field and determines the IMF 
described by means of Eqs. (IT3| [TBI [TBl - EUl) . at and the field polarity are used to deal with 
the drift velocity (as discussed in Sect. Hj), which modifies the overall convection velocity 
[Eq. (jl])]. Drift contribution is relevant during low solar activity - e.g., for at < 30° (Sect. [3]) 
- and decreas es with increasing solar acti vity, at values are obtained from Wilcox Solar 

and are c alculated using two different models 



Observatory (IHoeksema 



1995 



WSO 2010) 



called "R" and "L" . iFerreira fc Potrieted f l2003l . 120041 ) suggested that "R" model accounts 



for GCR observations during periods of increasing solar activity (for instance, 1987.4-1990.0 
and 1995.5-2000.0), while "L" model accounts for periods of decreasing solar activity (for 
instance, 1990.0-1995.5 and 2000.0-2010.0). The implementation of "R" and "L" models in 
the current code is further treated in Sects. I7.2H7.2.11 Finally, the latitudinal dependence 
(e.g., see Sect. [3]) of the solar wind [Eq. (ITTP ] depends (at low solar activity) on the values 
(averaged over 27 days) of SW speed and on the ecliptic at Earth orbit. 

The time spend by the SW to cover the distance from the outer corona up to the bound- 
ary of the effective heliosphere can be expressed in units of the time ne eded for a sidereal 
rotation on the equator of the Sun (about 2 5 days, e.g., see pag e 77 of lAschwandenl 12006 
and also iBrajsa et al.ll200ll : for a survey see iRuzdjak et al.ll2005l ). For instance depending 
on the wind speed, on the ecliptic the SW spends the corresponding amount of time needed 
to complete from 12 up to 20 sidereal solar rotations to reach the outer boundary. In the 
present code, the effective heliosphere (with a radius of 100 AU) was subdivided in 15 spheri- 
cal regions. In each region, the parameters (e.g., SW speed, Kq, B^, at, etc.) are determined 
at the time of the solar wind ejection. 
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The Monte Carlo Code HelMod 



It is worthwhile to note that Eqs. ([H H]) can be analytically solved onl y treating a simpli- 
fied transport of GCRs through the hel iosphere (e.g., see Sect. 12. H and also lGleeson fc Axford 
19681 : ICaballero-Lopez fc Moraalll2004j ). Complex configurations regarding the transport in- 
side the h eliosphere were proposed using numerical methods, like finite-difi^erence integra- 



tion (e.g.. iBurger fc Potgieteiill989l ). 



As ir nplemented in the H e lMod codqj v e rsion 1 .5, the curre nt approach i) follow s 
that from IVamada et aP fll998h: ICervasi et aP fll999h : IzhaneJ fll999h : lAlanko et all fcoosh : 



Pei et al.l ( l2010bl ): IStrauss et al.l ( 1201 ll ) and ii) exploits a Monte Carlo technique to deter- 
mine the number density U (Sect. [2]) using the set of the approximated stochastic differential 
equations (SDEs) treated in Appendix |X] for a 2-D approximation (radial distance and co- 
latitude). For a) an IMF described by the standard Parker field [Eq. ( !T5|) ] and b) both 
solar wind and drift velocity in the region of WHCS radially directed (e.g., = and 
"^^Hcs.r = "i^Hcs)) the SDEs approximated in terms of the increments Ar, Afi{9), AT and At 
[with fi{9) = cos(6')] are (see Appendix |A|) : 



Ar 
AM 



dr ^ 



At + (Vsw + ^^dr,r + Vncs) + COrV^K^^At 



(24) 



At, 



AT 



2 arciVs^T 



At 



(25) 
(26) 



[see Equations (2-4) of lBobik et al.ll2011dl . see also lPei et al.ll2010bl and references therein]. For 
the IMF treated in Sect. [3] [Eqs. (fT8H20|)] the SDEs [Eqs. fE25tiA27ll ] can be approximated 
with [Eqs. flA3TM33l) ]: 



Ar 



^2 Q^^ rrJ 



d 



dfi{e) 



— " ""'^'^ 



r01 



At 



(27) 



^In the 2D-HelMod code version 1.0, the standard Parker field without drifts was implemented; in version 
1.2, the dependence on the particle drift was added; finally, in version 1.4, the Parker magnetic field was 
modified in polar zones. 
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AT 



]_d_ 

^2 Qj, 



2 aj-ciVs^T 



At. 



■At 



(28) 
(29) 



As discussed by lPei et al.l (l2010b[ ) [see also ( IStrauss et al.ll201l[ )]. the vector q = (r, /x, T) 
represents a so-called pseudopartide (see Appendix . Equations ( 127H29|) allow one to sim- 
ulate the time evolution of pseudoparticles from the outer boundary down to the inner helio- 
sphere. As treated by Achterberg and Krulls (1992), the number density U - or equivalently 
the differential intensity J [Eq. ([3])] - can be obtained from the density of pseudoparticles by 
averaging over many realizations of the SDEs. 

The procedure to integrate the SDEs is the following: 1) events are isotropically gene- 
rated on the outer border of the effective heliosphere; 2) each event is integrated over the 
time evolution of a pseudopartide and is processed forward-in-time until it reaches either the 
outer (inner) border of the effective heliosphere located at 100 AU (ri,) or the pseudopartide 
energy becomes lower than a minimum threshold (which depends on the set of experimental 
data taken into consideration), then a new particle is generated; 3) when a pseudopartide 
reaches a particular region (for instance that corresponding to Earth position) its injection 
energy, statistical-weight, etc. are recorded; 4) finally, the number density U results from 
the normalized distribution function obtained using a procedure from iPei et al.l (j2010bl ) (see 
Section 4.3 in this article). The forward- in-time approach allows one to reproduce rigorously 
processes occurring inside the heliosphere. 

In the present code. At varies as r'^/K„, thus allowing an increase of the accuracy in 
the inner heliosphere, but keeping the appropriate precision up to regions close to the outer 
border of the effective heliosphere. Fu rthermore, this condition e nsures that the diffusion 
process is dominant (see Section 4.1 of iKruells fc Achterberg||l994j ). 



7. Results 

The current modulation code (Sect. |6]) provides a modulated differential intensity for 
protons using a local interstellar spectrum (LIS) of protons. In the following, we will discuss i) 
the LIS used (Sect. I7.ip . ii) the comparison of simulated (modulated) differential intensities 
with those obtained from the measurements of BESS, AMS and PAMELA spectrometers 
during the solar cycle 23 (Sects. 17.2^ I7.2.ip and iii) the dependence of present results on the 
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Fig. 3.— Left: spectral index (7) obtained (see text) for AMS-1998, BESS-1998, BESS-2002 
and PAMELA-2006/08. Right: normalization constant (Jq) for BESS-1997, AMS-1998, 
BESS-1998, BESS-1999, BESS-2002 and PAMELA-2006/08. The dotted lines represent 
the values of the spectral index (7bph) and normalization constant ( Jq.bhp) of the BPH-LIS, 
respectively; the continuous lines represent the error-weighted averages of spectral index 
(7wa) and normalization constant (Jo,wa)- 



treatment of the heliosphere extension (Sect. 17.30 . Furthermore (Sect. UAL the simulated 
fluxes obtained with HelMod code are compared with (a nd found to reproduce the feature s 
of) the experimental data from Ulysses fast scan in 1995 (jSimpson. Zhang and Bamdll996l ). 



7.1. Local Interstellar Spectrum 



Recently, iHerbst et al.l (120101 ) reviewed different proton LIS's published in the litera- 
ture and determined that - as it can be seen in Figure 2(b) in that article - these spectra 
agree well with each other for proton energies above 10 GeV. For this comparison, they used, 
among others, the LIS from Burger, Potgieter and Heber (2000) (BPH-LIS) in the form 
of the approximated analytical expression from Usoskin and collaborators (2005). Over the 
past years, M oskalenko, Strong and col l a.borators using GALPROP p r ovided a LIS f or pro - 
tons [e. g., see 



Moskalenko et al 



see also Langnerl (2004); Langner et al 



J2OO2I) : Istrong fc MoskalenkJ toO^ : IXrotta et al.l toill ) 



LIS above 1 GV [e.g., see iTrotta et al 



(|2003[) ]: the latest calculation agrees with the BPH- 
boil\ ]]. It has to be remarked that the GALPROP 
spectrum is constrained by a few measured quantities (for instance, the B/C and other iso- 
topes and/or nuclei ratios), some of them will be (accurately) re-determined in the coming 
years using data from PAMELA and AMS-02 missions. 
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In units of (srm^sGeV) ^ ( Burger et al.lboool . see also lUsoskin et al.l I2OO2I ) the BPH- 
LIS is expressed as: 



0,BHP 



^-7BPH for R>7, 



Jbhp(T) = < 



exp[9.472 - 1.999 Ini? - 0.6938 (Ini?)^ 
+0.2988 {InRf - 0.04714 (Ini?)^] , 



(30) 



for i? < 7, 



with 



where 



R = R{T) 



P{T) 



(e.g., see Equation (4.94) in lLeroy &: Rancoitall201ll ) is the proton rigidity in GV with i?r,p = 
mpC^, Trip is the rest mass of protons in GeV/c^, T is the kinetic energy of proton in GeV, 
e is the electron charge, c is the speed of light, 7bph = 2.78 is the spectral index, Jo,bhp = 
1.9 X 10'* (srm^sGeV) is a normalization constant and, finally, Pq = 1 GV. 

Above (10-20) GeV the differential proton intensities are slightly or marginally affected 
by modulation. The BPH-LIS [first line of Eq. (!30l) ] was compared to experimental spec- 
tra available in the literature and collected during the solar cycle 23. These observations 
also account for da ta in the energy range where modulation is relevant, e.g ., AMS-1998 



fjAguilar et al.ll2002[ ) BESS-1998 [with data only in the range (20 -117) GeVl flSanuki et al. 



20001), BESS-2002 flHaino et all l2004f ) and PAMELA-2006/08 flAdriani et al.l l2011ah. In 



Fig.pi the spectral indexes (7) of AMS-1998 and PAMELA-2006/08 are those from flAguilar et al 



2002 



Adriani et al.ll2011b[ ). respectively; while for BESS-1998 and BESS-2002 the spectral 



indexes were obtained from a fit to the published data of the differential proton intensities. It 
has to be noted that the rigidity independent part of the spectral index found by PAMELA- 
2006/08 is 7PAMELA = 2.790 ± 0.008(stat) ± O.OOl(syst); Adriani and collaborators (2011b) 
found that the spectral index depends on rigidity as expressed in Equation (19) therein 
with a maximum variation of the order of the previously quoted uncertainties in the rigidity 
range (30-200) GV . Furthermore, the spectral index (2.79 ± 0.08) found by Caprice-1994 
( iBoezio et al.lll999[ ) is in agreement with those found by the experiments discussed in this 
section, but the quoted errors are larger. 

The normalization constants Jn (FiS- E D a) depend o n the set of experim e ntal observa- 
tions , e.g., BESS- 1997 dshikaze et alJ l2007h. BESS -1998 Jshikaze et al .ll2007t ISanu ki et al. 



2000'), AMS-1998 keuilar et al.ll2002h. BES S -1999 fishikaze et al.lboOTL BESS-2002 Jflaino et al 
2004 ) and PAMELA-2006/08 jAdriani et al.l boilaf) and b) were obtained from a fit using 
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7bph as spectral index to the experimental data. For BESS-2000 (IShikaze et al.l 120071 ) . the 
experimental observations did not exceed the 21.5 GeV, i.e., an energy region of proton diffe- 
rential intensity which might (marginally) still be affected by modulation in a period of high 
solar activity; thu s, the normalizat ion constant used for these data was the one obtained 



from BESS-2002 (IHaino et al.l 120041 ) data. 



The weighted averages of both the spectral index (7) and normalization constan t (Jo) 



and th eir errors were determined following the procedure indicated at pages 14-15 of IPDB 



(120101 ). The error- weighted averages found are 

7wa = 2.783 ± 0.009 (31) 

and 

Awa = (1.76 ±0.01) X 10^ (srm^sGeV)"^ (32) 

7wa is well in agreement with that (7bph) suggested by Burger, Potgieter and Heber (2000) 
[Eq. (15U]) ]. Jo,wa and 7wa are represented with the continuous lines in Fig. [31 in the same 
figure the dotted lines refer to the values of the BPH-LIS [Eq. (15U]) ]. It has to be rernarked 



that the value of Jo found from a fit to Caprice-1994 data above 20 GeV (IBoezio et al.lll999l ) 



is 1.44 ± 0.02: this value differs by more than 5 standard deviations from Jo,wa [Eq- (l32|) ]. 

In Sects. I7.2H7.2.11 using the current modulation code the observed proton spectra are 
compared with the modulated differential intensities obtained from an interstellar differential 
(per unit of kinetic energy) proton intensity [JHeiMod(7')] given by 



JneiModiT) = Jbhp(T) ( ] [srm^sGeV] \ 

V-'0,BHP/ 



(33) 



JHeiModiT) keeps the same spectral index for P(T) > 7GV as in Eq. (1301) and linearly 
depends on Jq, which accounts for the slight absolute fluxes variation among observations. 



7.2. Comparison with Observations Obtained During Solar Cycle 23 

We used the present code for quantitative comparisons [using Eqs. ( IMl 1351) ] with expe- 
rimental data (discussed later in this section) collected during solar cycle 23, in periods with 
high solar activity, i.e., when the solar magnetic field becomes increasingly complex and less 
dipolar (Sects. [21 [3]). This code allowed us to investigate how the modulated (simulated) 
differential intensities are affected by the i) particle drift effect (Sects. [21 [Ij), ii) polar en- 
hancement of the diffusion tensor along the polar direction {K±g) [Eq. ([8])] and, finally, iii) 
the value of tilt angles (at) calculated following the approach due to "R" and "L" models 
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"L" model "R" model no drift diagonal approx. scalar approx. 



BESS-1999 


8.7 


8.0 


14.6 


32.0 


29.7 


BESS-2000 


16.2 


15.8 


13.0 


23.6 


26.7 


BESS-2002 


12.7 


15.0 


12.2 


34.8 


33.2 



Table 2: For BESS-1999, BESS-2000 and BESS-2002, r^RMs (in percentage) obtained from 
Eq. with enhancement of the diffusion tensor along the polar direction using "L" and 
"R" models for the tilt angle and for no drift approximation, diagonal approximation and, 
finally, scalar approximation (see text). 







"L" model 


"R" model 


no drift 


diagonal approx. 


scalar approx. 


BESS- 


-1999 


6.8 


8.1 


24.3 


30.5 


30.2 


BESS- 


-2000 


11.3 


10.2 


10.8 


26.2 


26.2 


BESS- 


-2002 


13.0 


15.7 


12.7 


33.9 


33.2 



Table 3: For BESS-1999, BESS-2000 and BESS-2002, //rms (in percentage) obtained from 
Eq. without any enhancement of the diffusion tensor along the polar direction using "L" 
and "R" models for the tilt angle and for no drift approximation, diagonal approximation 
and, finally, scalar approximation (see text). 



[Sect. Eland (IHoeksemal Il995l : IWSOll2010l )]. The magneti c field is modifi e d wit h respect to 



Parker's magnetic field in the polar region as proposed by lKota &: Jokipiil (119891 ) (Sect. |3]). 



The effects related to particle drift were investigated (a) via the suppression of the drift 
velocity - i.e., under the assumption that Ka = (Sect. H]), thus no drift convection was 
accounted for -, (b) in a pure diffusion approximation with a diagonal diffusion tensor (termed 
diagonal approximation), where K„ = JC and Kgg = p^JC (Sect. [2]) and, finally, (c) in a pure 
diffusion approximation with components both equal to /C (called scalar approximation) [as 
in Eq. (|9])]. The case (a) accounts the hypothesis that magnetic drift convection is almost 
completely suppressed during solar maxima. In addition, for cases (b) and (c) one allows to 
assume that the diffusion propagation is independent of magnetic structure. 

Each modulated (simulated) differential intensity was obtained using a diffusion tensor 
(Sects. |2| 12. U pland Appendix IXj) . whose elements depend on the actual value of the diffusion 
parameter Kq. Furthermore, the modulated spectra were derived from a LIS [Eqs. f pOj [33]) ] 
whose normalization constant ( Jq) depends on the experimental set of data (see discussion in 
Sect. I7.ip . In addition, these differential intensities were calculated 1) for a polar-increased 
value of K±g [Eq. ([8])] and also with K±g = K±r, and 2) accounting for particles inside two 
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Fig. 4. — Differential intensity determined with HelMod code (continuous line) compared to 
the experimental data of BESS-1999; the dashed line is the LIS (see text). 



heliospheric regions where solar latitudes are lower than |5.7°| and |30°|, respectively. As 
discussed in Sect. 12. in the present model Kp is assumed to be equal to the value of the 
rigidity jP) [Eg. (17|)] above proton kinetic energies of ~ 444 MeV (e. g., iGloeckler fc Jokipii 
19661 : loieeson fc Axfordl[l968l : IPerkJ[l987l : IPotgieter fc Le R,ouxlll994l ). However, it has to be 
remarked that a systematic investigation of its dependence below that value and the shape of 
low energy part of the LIS spectrum [Eqs. ( l30l 1331) ] was not attempted using the modulated 
intensities obtained from HelMod code. In fact, this investigation is likely to be carried out 
using the experimental data from accurate observations over a long duration, like those from 
the AMS-02 spectrometer which will allow one to reconstruct the particle trajectory. The 
reconstructed particle trajectory results in untangling GCRs coming from outside the magne- 
tosphere also at large geomagnetic latitudes (6m) where less energetic particles can enter the 
magnetosphere. For ins tance, inside highest geomagnetic reg ion with 1 < 6m < 1.1 radian 
[e.g., see Figure 2(c) in dAlcaraz et aPboOOah and Figure 8 in JSobik et al.ll2006h ] AMS-1998 
data indicate that i) the effective geomagnetic cut-off prevents primary protons (i.e., CR 
protons) from being fully observed with energies below ^ (0.5-0.6) GeV and ii) secondary 
particles largely contribute to the overall differential intensity. In addition, it has be noted 
that BESS observations were usually performed at large geomagnetic latitudes with 6m close 
to 1.13 radian. 



The past period of high solar activity was during the solar cycle 23; BESS collaboration 
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Fig. 5. — Differential intensity determined with HelMod code (continuous line) compared to 
the experimental data of BESS-2000; the dashed line is the LIS (see text). 



took data in the years 1999, 2000 and 2002 [see sets of data in (IShikaze et al.ll2007l )]. These 
data were compared with those obtained by means of HelMod code using the error-weighted 
root mean square (//rms) of the relative difference (rj) between experimental data (/exp) and 
those resulting from simulated differential intensities (/sim)- For each set of experimental 
data and above described approximations and/or models, we determined the quantity: 



with 

/sim(2~i) /exp(^j) /oir\ 

Vi = r~{T\ — ' ^^^^ 

Jexpy-l-i) 

where Tj is the average energy of the i-th energy bin of the differential intensity distribution 
and the errors including the experimental and Monte Carlo uncertainties; the latter 

account for the Poisson error of each energy bin. The simulated differential intensities are 
interpolated with a cubic spline function. 



In Tables [21 |3l the values of the parameter ?7rms (in percentage) are shown; they were 
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Fig. 6. — Differential intensity determined with HelMod code (continuous line) compared to 
the experimental data of BESS-2002; the dashed line is the LIS (see text). 



obtained in the energy rang^ f rom 444 Me V up to 30 GeV u sing "L" and "R" models for 
the tilt angle (at) [Sect. [5] and (jlloeksema 1995 : WSO 2010 )]. for no drift approximation, 
diagonal approximation and scalar approximation (approximations discussed previously in 
this section), finally with (Table [2]) and without (Table [3]) the enhancement of the diffusion 
tensor along the polar direction {Kj_g) [Eq. ([H])]. The simulated differential intensity were 
obtained for a heliospheric region where solar latitudes are lower than |30°|. From inspection 
of Tables [2] and [3l one can note that i) the no drift approximation is better appropriate 
than diagonal and scalar approximations, ii) the "L" model for calculating the values of tilt 
angle (at) is slightly to be preferred to "R" model (although the overall differences between 
these two models are marginal), iii) the results obtained accounting for drift effects using 
tilt angles from "L" model are better in agreement with experimental data with respect to 
the no drift approximation and, finally, iv) the minimum difference with the experimental 
data occurs when K±g = K±r is assumed independently of the latitude (Table [3l see first 
column of the left-hand side). In addition, the results obtained for a heliospheric region where 
solar latitudes are lower than |5.7°| exhibit a behavior similar to those lower than |30°|, but 
with values of t/rms (in percentage) larger by about several percents. In Figs. HI [5l [6] the 



^Above 30 GeV the differential intensity is marginally (if at all) affected by modulation. 
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"L" model 


"R" model 


no drift 


diagonal approx. 


scalar approx. 


BESS-1997 9.2 


17.7 


10.4 


9.5 


17.6 


AMS-1998 4.6 


7.9 


12.9 


5.4 


17.3 


BESS-1998 9.1 


14.1 


9.3 


4.7 


13.6 


PAMELA-2006/08 7.1 


13.4 


5.9 


17.5 


52.5 


Table 4: For BESS-1997, AMS-1998, BESS-1998 and PAMELA-2006/08, i^rms (in per- 
centage) obtained from Eq. (l34l) with enhancement of the diffusion tensor along the polar 
direction using "L" and "R" models for the tilt angle and for no drift approximation, diagonal 
approximation and, finally, scalar approximation (see Sect. 17. 2|). 


"L" model 


"R" model 


no drift 


diagonal approx. 


scalar approx. 


BESS-1997 13.4 


20.6 


14.2 


11.13 


12.0 


AMS-1998 6.1 


11.3 


11.4 


6.0 


3.7 


BESS-1998 11.1 


17.7 


7.3 


4.1 


7.1 


PAMELA-2006/08 11.0 


24.7 


5.4 


12.3 


30.6 



Table 5: For BESS-1997, AMS-1998, BESS-1998 and PAMELA-2006/08, t^rms (in per- 
centage) obtained from Eq. f lM|) without any enhancement of the diffusion tensor along the 
polar direction using "L" and "R" models for the tilt angle and for no drift approximation, 
diagonal approximation and, finally, scalar approximation (see Sect. I7.2p . 



differential intensities determined with HelMod code are shown and compared with the 
experimental data of BESS-1999, BESS-2000 and BESS-2002, respectively; in the same 
figures, the dashed line is the LIS [Eqs. f pOj [33]) ] with normalization constants Jq treated in 
Sect EH These modulated intensities are the ones calculated for a heliospheric region where 
solar latitudes are lower than |30°|, using K±0 = K±r independently of the latitude and 
including particle drift effects with the values of tilt angle from the "L" model. 

Finally, it has be concluded that the present code combining diffusion and drift mecha- 
nisms is suited to describe the modulation effect in periods with high solar activity (e.g., see 
Ferreira fc Potgieterll2004] : iNdiitwani et al.l 120051 ). 



7.2.1. Periods not Dominated by High Solar Activity 



In periods where the solar activity is no longer at maximum, the solar magnetic field 
becomes increasingly dipolar (Sects. [21 13]). We used the present code to compare the simu- 
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Fig. 7. — Differential intensity determined with HelMod code (continuous line) compared to 
the experimental data of BESS-1997; the dashed line is the LIS (see text). 



lated differential intensities with experimental data obtaine d during periods no t dominated 
by high solar activity in the solar c ycle 23, i.e., BE SS-1 997 Jshikaze et al.lboOTh . AMS-1998 
jAguilaret al.l 120021) . B ESS-1998 Jshikaze et al.l I2OO7I : ISanuki et al.lboool ) and PAMELA- 



2006/08 (lAdriani et al.l l2011al ). As discussed in Sect. \7.2\ the simulated spectra were cal- 
culated including the effects due to particle drift - expected to be relevant (Sects. El H]) - 
with the value of ti lt angles (at.) calculated fol lowing the approach due to "R" and "L" 
models [Sect.[5]and ( lHoeksemalll995l : IWSOll2010l )]. with and without the polar enhancement 
of the diffusion tensor along the polar direction {K±0) [Eq. (|8])]. Similarly to the treatment 
for periods with high solar activity (Sect. 17. 2p . the effects related to particle drift were also 
investigated (a) via the suppression of the drift velocity (no drift), (b) with the diagonal 
approximation and, finally, (c) with the scalar approximation. 

In Tables m and |5], the values of the parameter t^rms (in percentage) are shown. They were 
obtained in the energy ra nge from 444 MeV up to 30 G eV using "L" and "R" models for the 
tilt angle (at) [Sect. Eland ( lHoeksema]ll995l : IWSOll2010l )] . for no drift approximation, diagonal 
approximation and scalar approximation (approximations discussed in this Sect l7.2l) . finally 
with (Table H]) and without (Table [5]) the enhancement of the diffusion tensor along the 
polar direction {Kj_g) [Eq. ([8])]. The simulated differential intensity were obtained for a 
heliospheric region where solar latitudes are lower than |5.7°|. From inspection of Tables H] 
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Fig. 8. — Differential intensity determined with HelMod code (continuous line) compared to 
the experimental data of AMS-1998; the dashed line is the LIS (see text). 

and [5], one can note that i) the diagonal approximation is better appropriate than no drift 
and scalar approximations, ii) the "L" model for tilt angles (at) is slightly to be preferred 
to "R" model and, finally, iii) the minimum difference with the experimental data occurs 
when the enhancement of the diffusion tensor along the polar direction (K^g) [Eq. ([8])] is 
taken into account (Table HI see first column of the left-hand side). In addition, the results 
obtained for a heliospheric region where solar latitudes are lower than |30°| exhibit a behavior 
similar to those lower than |5.7°|, but with values of r/RMS (in percentage) larger by about 
several percents. In Figs. [71-[T0| the differential intensities determined with HelMod code are 
shown and compared to the experimental data of BESS-1997, AMS-1998, BESS-1998 and 
PAMELA-2006/08, respectively; in the same figures, the dashed line is the LIS [Eqs. (I5UI I55]) ] 
with normalization constants Jq treated in Sect l7.1[ These modulated intensities are the ones 
calculated for a heliospheric region where solar latitudes are lower than |5.7°|, using the 
enhancement of the diffusion tensor along the polar direction (K^q) [Eq. ([8])] and including 
particle drift effects with the values of tilt angle from "L" model. 

Finally, it has be concluded that the present code combining diffusion and drift mecha- 
nisms is also suited to describe the modulation effect in periods when the solar activity is 
no longer at the maximum. 
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Fig. 9. — Differential intensity determined with HelMod code (continuous line) compared to 
the experimental data of BESS-1998; the dashed line is the LIS (see text). 



In Sects. 17.21 and 17.2.11 the simulated differential intensities were obtained from a LIS 
[described by Eqs. ( 15U| 132])] propagating through a spherical heliosphere with a radius of 
100 AU down to Earth. However, the physical dimensions of the heliosphere also depends on 
the speed of solar wind. In HelMod code, the simulated modulated intensities are determined 
by the properties of the diffusion tensor (Sects. O 12.11 |4]and Appendix |A]), whose elements 
are related to the actual value of the diffusion parameter. Kq acts as a scaling factor for the 
overall modulation effect. It was indirectly determined from neutron monitor measurements, 
thus, it is expected to be sensitive to the overall modulation effect (from the heliosphere 
boundary down to Earth), but almost independent of the variation of heliosphere dimensions. 

The radial distance of the heliosphere was varied from 80 up 120 AU. The corresponding 
simulated differential intensities were compared to the experimental data from BESS-2002 
[data collected during high solar activity (Sect. O])] and PAMELA-2006/08 [data collected 
when the solar activity was no longer large (Sect. I7.2.ip ]. i.e., when heliosphere is expected 
to be smaller or larger (and possibly no longer spherical) than 100 AU, respectively. 

The values of ?7rms [Eq- (13^ ] in percentage calculated for spherical heliospheres with 



7.3. Dependence on the Extension of Heliosphere 
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Fig. 10. — Differential intensity determined with HelMod code (continuous line) compared 
to the experimental data of PAMELA-2006/08; the dashed line is the LIS (see text). 



radii of 80, 90, 110 and 120 AU are shown in Tableland compared with those calculated 
with a radius of 100 AU (see Tables [3l H]). For BESS-2002, the sim ulated intensities were 
obtained i) using the "L" model for the tilt angle (at) [Sect. [5] and (IHoeksemal Il995l : IWSO 
20101 )]. ii) with Kx^o = K±r independently of the latitude and iii) inside a heliospheric region 



where solar latitudes are lower than |30°|. For PAMELA-2006/08, the sim ulated intensities 
were obtain ed a) using the "L" model for the tilt angle (at) [Sect. EJand (IHoeksemal Il995 



WSOll2010l )]. b) with an enhancement of the diffusion tensor along the polar direction {K^g) 
[Eq. dH])] and c) inside a heliospheric region where solar latitudes are lower than |5.7°|. From 
inspection of Table El one can remark that, within 2.3%, the simulated differential intensities 
for spherical heliospheres with radii of 80, 90 and 110 AU are compatible with that with a 
radius of 100 AU; slightly larger values of //rms were obtained for a spherical heliosphere with 
a radius of 120 AU. 

The sensitivity of this approach was estimated from the differences of the simulated 
intensities with radii of 80, 90, 110 and 120 AU with that with a radius of 100 AU for protons 
with energies above 30 GeV, i.e., for an energy region in which the spectrum is unaffected by 
modulation and, thus, no difference is expected. For this purpose, we defined the quantity 
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[see also Eqs. f p^ [35 



'7RMs,h - y — ^ ^^^2 — (36) 

with 



_ /hl^i) — flOOAv{Ti) , . 

JlOOAuUij 

where fh{Ti) is the differential intensity of z-th energy bin (above 30 GeV), crj^,i,h is the error 
due to Monte Carlo uncertainties for i-th energy bin, /looAu(^j) is the differential intensity 
computed with a radius of 100 AU and, finally "h" indicates 80, 90, 110 and 120 AU. //RMS.h 
resulted equal to about 2.3% for heliospheres with 80, 90, 110 and 120 AU. Thus, the mo- 
dulated intensities for heliospheres with radii of 80, 90, 100, 110 AU (and also 120 AU for 
BESS-2002) are in agreement among them and experimental data within the present sensi- 
tivity of about 2.3% of the current approach; at 120 AU the simulated intensity is marginally 
non compatible with that obtained with 100 AU for PAMELA-2006/08. These results indi- 
cated that, as expected, the diffusion parameter almost accounts for effects related to the 
variation of the physical dimensions of the heliosphere within the present approximations. 



7.4. Dependence on Heliospheric Latitude 



Observations made by the Ulysses spacecraft ( ISimpson et al.l 119921 ) in the inner helio- 
sphere could determine a latitudinal dependence of GCR (mostly protons) intensity with an 
equatorial southward offset minimum and a North polar excess. This dependence was dis- 
cussed also in terms of modulatio n models which were including particle drift effe cts (e.g., see 



Simpsonlll996l : lHeber et al.lll998l ). For protons with energies larger than 100 MeV. ISimpson. Zhang and Bam 
fll996[ ) expressed their results in terms of the solar latitude and found that i) the latitudinal 
gradient is ~ (0.33 ±0.02)% deg~^, ii) the c ounting rate minimum is ne a rly co nstant in a lat- 
itudinal r egion of y — (15°-5°) [Figure 2 o f ISimpson. Zhang and Bamd (Il996[ )] at ~ 1.35 AU 
[e.g., see (ISimpsonlll996t iHeber et al.lll998l )] and iii) the rate at the minimum is about ~ 80% 



80 AU 90 AU 100 AU 110 AU 120 AU 



BESS-2002 


14.5 


12.00 


12.2 


13.0 


14.5 


PAMELA-2006/08 


5.7 


6.1 


7.1 


7.7 


10.8 



Table 6: ?7rms (in percentage) calculated for a spherical heliosphere with radius of 80, 90, 
110 and 120 AU and compared with those obtained with a radius of 100 AU: for BESS-2002 
(see Table ED, for PAMELA-2006/08 (see Table H. 
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Fig. 11. — i?iat calculated at 1 AU as a function function of solar latitude for a) Te > 444 MeV, 
b) Te > 600 MeV, c) > 1200 MeV and d) > 2100 MeV: left-hand side during 1995 
(A > 0), right-hand during 2007 (A < 0). 



with respect to that at ~ 80°. iFerreira et al.l (120031 ) have shown that the latitudinal dip - 
found with the Ulysses fast scan - can be reproduced in a model using the Parker standard 
field and a polar enhancement of the diffusion tensor. 

Using HelMod code, we could investigate the latitudinal dependence of the differential 
intensity at 1 AU, above 444 MeV (as so far treated) up to 200 GeV. The heliosphere was 
subdivided in 20 regions equally spaced with respect to the co-latitudinal parameter fj,{9) 
[Sect.E]. The total fluxes obtained in each region were divided by the maximum flux occurring 
at the North pole, thus, Ri^t represents the normalized flux (to that at the North pole) as 
a function of the co-latitude. In addition, the values of i?iat were calculated for periods of 
opposite magnetic polarities and compatible with Ulysses pole-to-pole fast scans, i.e., for the 
years 1995 with A > and 2007 with A < 0. Ri^t can be equival ently expressed as a function 
of the solar latitude for a comparison with the results obtained by lSimpson. Zhang and Bame 
(119961 ). i?iat as a function the solar latitude is shown in Fig. [11] for the year 1995 (left- 
hand side) and 2007 (right-hand side). i?iat is also shown as a function of the minimum 
kinetic energy accounted for protons (Te), i.e., a) Tg > 444 MeV, b) Tg > 600 MeV, c) 
Te > 1100 MeV and d) Tg > 2100 MeV. By inspection of Fig. [IH for the year 1995 one 
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can remark that i?iat has 1) a latitudinal gradient of (0.23 ± 0.01)% deg~^, 2) an equatorial 
southward offset minimum in the latitudinal region ^ — (18°-5°) with Tg > 444 MeV and 3) 
at this minimum, the flux is about ~ 80% of that at the North pole. Thus, the simulated 
fluxes re produce the features of the experim ental dat a (above 100 MeV) exhibited in Figures 2 
and 3 of lSimpson. Zhang and Bamd (119961 ) [see also (IHeber et al.ll2008l )] regarding the period 
of 1995 Ulysses fast scan. However with increasing Tg, the latitudinal gradient decreases 
(Fig. [m left-hand side). In the year 2007, a similar minimum is exhibited for Tg > 444 MeV 
with a ~ 10% North-South poles asymmetry; with increasing Tg, this asymmetry gradually 
disappears and the flux reduction on the equatorial region is less pronounced down to ~ 88% 
with Tg > 2100 MeV (Fig. HH right-hand side). 

It is worthwhile to note that in HelMod code the magnetic field structure is treated 
similarly in North and South hemis phere approximating Parker's magnetic-field with that 
suggested by iKota fc Jokipiil ( 1l989l ) (Sect. |3]). However, the current 2-D model uses the 
complete 2x2 diffusion tensor (see Sects. [3l [6] and Appendix [X]) which contains both sym- 
metric and antisymmetric components in the off-diagonal terms. The symmetric component 
of the off-diagonal terms [Eq. (IA13P ] is determined by the divergenc e-free IMF used , which 
exhibits a latitudinal component arising from the modification by iKota fc Jokipiil (119891 ) 
[Eqs. (|T5| [T8| [T9|) ]. In the framework of the present 2-D model, the North polar excess and 
equatorial southward offset minimum shown in Fig. [11] originate from the non-zero symme- 
tric component of the off-diagonal terms. The actual extension of the dip is related to both 
the enhancement of the diffusion tensor in the polar regions [Eq. ([H])] and drift effects. 



8. Conclusions 

A systematic investigation of the solar modulation effect on the propagation of cosmic 
protons through the heliosphere down to the Earth was carried out comparing experimental 
observations performed during the solar cycle 23 and simulated differential intensities ob- 
tained using HelMod code. The simulated spectra were derived from a LIS [Eqs. (ISUj ES])], 
whose normalization constant ( Jq) depends on the experimental set of data (see discussion in 
Sect J7.ip . The stochastic 2D Monte Carlo (HelMod) code includes i) a fully treated diffusion 
tensor with symmetric and antisymmetric off-diagonal elements, b) a diffusion parameter 
which is a function of the intensity of solar activity and varies with solar polarity and phase 
(Sects. 12. H and c) a magnetic-field A vhich is modified with respect to Parker's magnetic 
field in the polar region as proposed by iKota Sz Jokipii Jl989h (Sect. ED. 



For observations performed during high solar activity, the simulated intensities (found 
with a better agreement to experimental data) were obtained i) using the "L" model for the 
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tilt angle (at) [Sect.[5]and flHoeksemalllQQSl : IWSOll2010h ]. ii) with K^e = K^r independently 
of the latitude and iii) inside a heliospheric region where solar latitudes are lower than 
|30°|. For observations performed when solar activity is no longer at the maximum, the 
simulated intensities (found with a better agreement to e xperimental data) were obta ined 
a) using the "L" model for the tilt angle (at) [Sect. [5] and JHoeksemallioOsi IwSOlboioh ]. b) 
with an enhancement of the diffusion tensor along the polar direction {K±e) [Eq. ([8])] and c) 
inside a heliospheric region where solar latitudes are lower than |5.7°|. 

In addition (within 2.3%), the simulated differential intensities for spherical heliospheres 
with radii of 80, 90 and 110 AU (and also 120 AU for BESS-2002) are compatible with that 
with a radius of 100 AU; a slightly lower agreement was obtained for a spherical heliosphere 
with a radius of 120 AU for PAMELA-2006/2008. These results indicated that, within the 
present approximations, the diffusion parameter almost accounts for effects related to the 
variation of the physical dimensions of the heliosphere. 

The simulated modulated spectrum determined for the year 1995 exhibits a latitudinal 
gradient of (0.23 ±0.01)% deg~^, an equatorial southward offset minimum in the latitudinal 
region — (18°-5°) with T^. > 444 MeV and at this minimum the flux is about ~ 80% of 
that at North pole. Thus, the simul ated fluxes reproduce the feature s of the experimental 
data from Ulysses fast scan in 1995 ( ISimpson. Zhang and Bamelll996l ). 



Although the treatment is highly simplified with respect to the complexity of physical 
mechanisms responsible for modulation effects, the overall satisfactory agreement found al- 
lows one to remark that the choice of parameters regarding the structure of IMF, diffusion 
tensor, diffusion parameter and tilt angle is almost appropriate to describe the experimental 
data. Finally, the experimental data from accurate observations over a long duration (like 
those from the AMS-02 spectrometer) will allow one to undertake a deeper systematic inves- 
tigation of solar modulation effects over a period longer than a solar cycle. Thus, possibly, 
further advancements can be put forward in the present approximations on the transport of 
GCR's through the heliosphere, for instance those at low rigidities, the spatial and rigidity 
properties of diffusion tensor. 
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Appendix 



A. Diffusion Tensor and Stochastic Differential Equations 

In a reference frame witli tlie 3rd coordinates along the avera ge magnetic field, the 
matrix of the diffusion tensor used in Eqs. ([H Sj) is given by (e.g., see |Jokipiilll97l[ ): 



ik 



Ka 




-Ka 

K^e 







Kn 



In heliocentric spherical coordinates (r, 6, (p), for instance those used in Eqs. flA2511X27l). the 



matrix elements of the 3x3 tensor are found, for instance, in Equation (17) of iBurger et al. 



( 120081 ). In a 2-D approximation, the matrix elements of the resulting tensor are: 



Krr = K±e sin^ ^ + cos^ ^ [K\\ cos^ ip + sin^ ip) , 

Kgg = K±g cos^ ^ + siu^ ^ (^K\\ cos^ + K^^ sin^ ip) ) 

Kr-e = —Ka sin + sin^ cos ^ | cos^ -ip + K±r sin^ -ip — K±q^ 

Kgr = Ka sin ip + sin^ cos ^ (^K\ \ cos^ ip + K±r sin^ ip — K±g^ 



(Al) 
(A2) 
(A3) 
(A4) 



with tan?/; = -B^/ B^. + Bj) and tan^ = Be /Br [see Figure 6 of Burger et al.l (12008f l]. 
where ip is the spiral angle (for a standard Parker IMF isjiip reduces to Eq. (fT6|l and 
tan^ = 0), and with 



K\\=p3h{r,t)Kp{P,t) 

K±r = PkK\\, 
K±g = L{9)pkK\, 



(A5) 

(A6) 
(A7) 



where pk = 0.05 and l{6) is a step function that is 1 in equatorial region and 10 in polar 
region (Sect. [2]). The diffusion parameter fci(r, t) is replaced by Kq for an effective heliosphere 
of 100 AU (Sect. 12.11) . Furthermore, the matrix elements of the later tensor consist of a 
symmetric {Kf/.) and antisymmetric {Kpj.) part: 



K 



ik 



Kfk + K 



with the antisymmetric part related to drift velocity [Eq. (l22i) ]. v^, and treated in Sect. HI Fi- 
nally, using Eqs. (1 A 511 ASP . Eqs. (1A1HA4I) can be re- written as: 



Krr = K^r = /3KoKp{P,t) 



3B 



[L{e)pk sin2^ + cos^^(cos^?/' + pfc sin^^)] , (A9) 
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Kre 



r0> 
5 

re 



3B 



[L{e)pk cos^ ^ + sin^ ^( cos^ + pk sin^ tp)] , (AlO) 

(All) 
(A12) 



with 



f3KoKp{P,t) 



re 

A _ T^A 



3B 



{sin^ cos^ [cos^ ip + Pk sin^ ip ~ '•(^)Pfc] } ? (^13) 

(AM) 



Equations ([T], Hj) can be re-expressed in heliocentric spherical coordinates as: 



dU 
'dt 



r sm & ou \ or 



'90 



'09 



d 



r sin 9 d(j) 
1 dr'^VrU 



1 



^ ar r o9 r sm t/ 



asin^l/et/ 



r sin 6' 



(9^ 



r sin 9 



1 dr'^VdrU 1 dsm9vdeU 



1 <9wd,0f^ 



dr 
1 dr^Vr 

y2 Qj. 



r sin 9 



39 



1 dsm9Ve 
r sin 9 89 



rsm9 
1 dVs 



rsin^ d(j) J dT 



d 



— {arciTU), 



(A15) 



where U is the number density of GCRs (Sect. [2]) and T is the kinetic energy. In turn, in a 
2-D (radial distance and co-latitude) approximation, Eq. flAlSP can be re-expressed as: 



dU 



1 d 



d 



d 



1 dr'^VrU 1 



d 



d 



sin 9 



"^^09" J ' rsin^S^ 
d sin 9VeU 1 dr^VdrU 



sm9Kl—U + K^a^U 

or r 



d_ 
'89' 



1 d sin 9vdeU 



+- 



dr 

1 dr^Vr 

y.2 Qj. 



r sin^^ 



89 



dr 



r sin 9 



89 



1 dsin9Vo\ d 



dT 



r sin 9 89 



(A16) 



Let us define the variable p = cos{9), then we obtain 

09 = -(1 - pY'^-^dp. 



(A17) 
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In addition, one can introduce the function 

F = Ur^. 

Using Eqs. (IA16I - IAT8|) . for a SW radially propagating (i.e. K^.r = Kw) we find: 

_d_ 

dr 

_d_ 
_d_ 



(A18) 



dF 
'dt 





d 


-1; 






dr 







F d 



J.2 Qj, 



d_ 

' df 



dr 



d_ 

dr 

2\ 



3r2 



2drdr^ ^ 2drdfi 



^ld_d_ 
2 djj, dr 



2 dfi dfi 



(A19) 



Furthermore, following the treatment discussed in Sections 4.3-4.3.5 of (IGardinerl Il985h . 
one can a) express the Fokker-Plank equation involving F - which, in turn, is a function of 
q = (r, /X, T) - as: 

d ^ \ ^ d , , , , 1 \ ^ d d 



{[D(q,t)],,F} 



(A20) 



with D = LL^ and b) obtain the equivalent set of differential equations 

dq = A(q, t)dt + L(q, t) dW(t), 



(A21) 



wher e A(q, t)dt accounts for the so-called advective processes (e.g., iKruells fc Achterberg 



19941 ). L(q, t) dW(t) is the stochastic term co ntaining dW(t) which is the increment of the 



so-called Wiener process (e.g.. Section 4.3 of IGardinerl Il985l ) . Equations ( IA21|) are termed 
stochastic differential equations (SDEs). 

Furthermore, one can note that i) the first right-hand term of Eq. (IA20I) is equal to 
those included in the first three lines of Eq. (IA19I) and ii) the second right-hand term of 
Eq. f lA20p is equal to those included in the fourth and fifth line of Eq. flA19p . Thus, using 
Eqs. flAT9l IA20|) one derives: 



Sr^ dr 



■"d,e 



(A22) 
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and 



D 



r 



(A23) 



As discussed by lPei et al.l (l2010bl ) - see Section 2 therein -, the matrix D can be downgraded 
to a two-by-two matrix, because second order acceleration mechanisms are not considered 
in Eqs. dHll]). 



As already shown by lPei et al.l (j2010bl ) - see Appendix B therein -, the matrix L is not 
unique. However, there is only one positive definite, i.e.. 



1/2 



2/^1, (1-^2) 



1/2 



1/2 



(A24) 



Finally, for a 2-dimensional model (like that treated in Sects. |3l H]), from Eqs. (]A21t IA22P 
and using Eq. (]A24I) one finds the following SDEs: 



dr 



r2 Or ^ dfi 



dt + 



0.5K 



s 



(A25) 



d/i 



1| (rKt^) + I 



dt + 



dT 



2 a^c\Vsy,T 



(A26) 
(A27) 



with dWi [i = r, fi{6)] the increment of the Wiener process. It has be remarked that the 
above usage of Eq. (1A24I) ensures imaginary terms do not appear in Eqs. ( ]A25tiA30l) . For 
an IMF described by a standard Parker field [Eq. (1151) ] requiring K^g = 0, the above SDEs 
reduce to 



dr 
d/i 



d_ 

dr 



dt + (V;„ + VA,r) dt + ^/2KS^ dWr 



(A28) 
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dT 



2 arelVswT 



dt 



(A29) 
(A30) 



Gardiner (1985) - see Section 4.3.1 therein - demonstrated that, following a Euler- 
Cauchy procedure, the SDEs can be approximated in terms of the increments Ar, Afi and AT 
occurring after a time At has elapsed. The corresponding increment of the Wiener Process is 
given by cUjA/At [with i = r, /i(^)], where ooi is a random number follo wing a Gaussian distri- 
butio n with a mean of zero and a stand ard deviation of one (e.g., see iKruells fc Achterberg 



1994 and appendix A of lPei et al.ll2010bf l. As a consequence, the SDEs [Eqs. flA25tiA27p 
be approximated by: 



can 



Ar 



At 



At. 



Afi 



d_ 



v± 



'2^-1,(1-^2) 
+0Ju.\l — . ^ ' At. 



AT 



2 a^ciVs^T 



At. 



(A31) 



At + 



(A32) 
(A33) 
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